Read in Data for Myeloarchitecture Characterization

Final study sample

participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv") #215 scans from 140 subjects. csv generated by /sample_construction/finalsample_7Tmyelin.Rmd

Glasser regions list

glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.plotting <- glasser.regions
glasser.plotting$cortex <- "cortex"
glasser.plotting$cortex[glasser.plotting$orig_parcelname %in% glasser.snr.exclude$orig_parcelname] <- "exclude"
glasser.plotting <- glasser.plotting %>% select(label, cortex)
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "rh_???", cortex = "exclude"))
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "lh_???", cortex = "exclude"))
glasser.plotting <- glasser.plotting %>% filter(label != "lh_L_10pp")

Depth-dependent R1 in HCP-MMP regions (Developmental study sample)

#R1 in all HCP-MMP regions at all 11 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.R

#R1 by region matrices in superficial, middle, and deep cortex
myelin.compartments.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/compartmentsR1_glasseratlas_finalsample.RDS")

Depth-dependent R1 in HCP-MMP regions (MICA-PNI high resolution sample)

#R1 by region matrices at all 14 cortical depths
myelin.glasser.micapni <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/MICA-PNI/depthR1_glasseratlas_MICAPNI.RDS") #generated by /surface_metrics/MICA-PNI/extract_depthdependent_R1.R

Myelin maps

#T1w/T2w Ratio, from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
T1wT2w.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/T1T2ratio/T1T2ratio_MeanMap_Baum2022_allages.pscalar.nii")
T1wT2w <- data.frame(T1wT2w = T1wT2w.cifti$data, orig_parcelname = names(T1wT2w.cifti$Parcel))

#Magnetization Transfer, from Hettwer et al., 2024, Nat Comm https://www.nature.com/articles/s41467-024-50292-2
MTsat <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/NSPN_MTsat/NSPN_baseline_average_MT.csv") %>% t() %>% as.data.frame() %>% purrr::set_names("MTsat")
MTsat$orig_parcelname <- paste0(rownames(MTsat), "_ROI")
MTsat$orig_parcelname <- gsub("\\.", "-", MTsat$orig_parcelname)

#Myelin Basic Protein Gene Expression, from Wagstyl et al., 2024, eLife https://elifesciences.org/reviewed-preprints/86933
MBP.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/AHBA_magicc/expression_maps/source-magicc_desc-MBPexpression_space-fsLR_den-32k.pscalar.nii")
MBP.lh <- data.frame(MBP = MBP.cifti$data, orig_parcelname = names(MBP.cifti$Parcel))
MBP.rh <- MBP.lh #reflect gene expression data across hemispheres for a whole brain map
MBP.rh <- MBP.rh %>% mutate(orig_parcelname = paste0("R_", substring(orig_parcelname, 3)))
MBP <- rbind(MBP.rh, MBP.lh)

Studying Cortical Depths Dominated by Signal from Cortical Gray Matter

#Compute average frontal cortex volume fraction at each depth for final study participants
cortexpve <- read_parquet("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/7T_BrainMechanisms_cortexpve.parquet") #parquet with cortex fraction measures
cortexpve <- cortexpve %>% filter(atlas == "glasser") #just for glasser atlas 
cortexpve <- cortexpve[!(cortexpve$StructName %in% glasser.snr.exclude$orig_parcelname),] #just for good regions
cortexpve <- cortexpve[cortexpve$StructName %in% glasser.frontal$orig_parcelname,] #just for frontal lobe

extract_depth_cortexpve <- function(mymeasure){
  measure.df <- cortexpve %>% select(subject_id, session_id, StructName, all_of(mymeasure))
  names(measure.df) <- c("subject_id", "session_id", "orig_parcelname", "cortex.fraction")
  depth.cortexpve <- merge(measure.df, participants, by = c("subject_id", "session_id")) %>% select(cortex.fraction) %>% colMeans(.$cortex.fraction) #get cortex.fraction for final study participants only, and then average for this depth
  return(depth.cortexpve)}

pve.measures <- list("Mean_cortexpve.1.00%", "Mean_cortexpve.0.9%", "Mean_cortexpve.0.8%","Mean_cortexpve.0.7%","Mean_cortexpve.0.6%","Mean_cortexpve.0.5%","Mean_cortexpve.0.4%","Mean_cortexpve.0.3%","Mean_cortexpve.0.2%","Mean_cortexpve.0.1%", "Mean_cortexpve.0.0%") #measures to extract data for 
cortexpve.frontallobe <- lapply(pve.measures, function(x) { 
  extract_depth_cortexpve(x)}) 
cortexpve.frontallobe <- do.call(rbind, cortexpve.frontallobe) %>% as.data.frame()
cortexpve.frontallobe$depth <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k")
cortexpve.frontallobe$depth <- factor(cortexpve.frontallobe$depth, ordered = T, levels = c("k", "j", "i", "h", "g", "f", "e", "d", "c", "b", "a"))
kable(cortexpve.frontallobe)
cortex.fraction depth
0.7103157 a
0.8388772 b
0.9281964 c
0.9749537 d
0.9922900 e
0.9960812 f
0.9925421 g
0.9710414 h
0.9071242 i
0.7684002 j
0.5157232 k
cortexpve.frontallobe %>% 
  ggplot(aes(x = cortex.fraction, y = depth, color = depth)) +
  geom_point(size = 2) +
  theme_classic() +
  xlab("\nCortex fraction") +
  ylab("Cortical depth\n") +
  scale_color_manual(values = c("#b7b7b7", "#b7b7b7", "#004a38", "#006458", "#0093b0", "#8aabd5", "#c0bfe3", "#e0d5ec", "#f2e7f4", "#b7b7b7", "#b7b7b7")) +
  theme(
  legend.position = "none",
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Depthplot_cortexfraction_frontal.pdf", device = "pdf", dpi = 500, width = 2, height = 1.6)
#Only study 7 depths where cortex fraction > 90%
myelin.glasser.7T <- myelin.glasser.7T[3:9]
names(myelin.glasser.7T) <- c(1, 2, 3, 4, 5, 6, 7) #1 is superficial (20% thickness); 7 is deep (80% thickness)

R1 Indexes Cortical Myelin Content Across Cortical Regions

myelin.colorbar <- c("#edf6ff", "#c2d4ed", "#809dc4", "#345c91", "#0e4669")

Whole-brain map of R1

Regional average R1 in HCP-MMP regions

#Calculate across-depths average regional R1
R1 <- do.call(rbind, myelin.glasser.7T) #R1 at analyzed depths, 20%-80% of thickness
R1 <- R1 %>% select(contains("ROI"))
R1 <- colMeans(R1) %>% as.data.frame()  %>% set_names("mean.R1")
R1 <- R1 %>% mutate(orig_parcelname = row.names(R1))

#Save out lh and rh mean R1 data 
R1.lh <- R1$mean.R1[1:180]
write.table(R1.lh, "/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/R1.lh.glasser.txt", col.names = F, row.names = F, quote = F)
R1.rh <- R1$mean.R1[181:360]
write.table(R1.rh, "/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/R1.rh.glasser.txt", col.names = F, row.names = F, quote = F)
df.list <- list(R1, T1wT2w, MTsat, MBP) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.9, height = 2.7)

Myelin map PCA

#run PCA
maps.PCA <- myelin.maps %>% select(label, T1wT2w, MTsat, MBP)
myelin.map.PCA <- maps.PCA %>% select(-label) %>% prcomp(., scale = T, center = T)
#variance explained
summary(myelin.map.PCA)
## Importance of components:
##                           PC1    PC2     PC3
## Standard deviation     1.5097 0.6641 0.52895
## Proportion of Variance 0.7598 0.1470 0.09326
## Cumulative Proportion  0.7598 0.9067 1.00000
#PC1
myelin.map.PC1 <- myelin.map.PCA$x[,1] %>% as.data.frame %>% purrr::set_names("PC1")
myelin.map.PC1$label <- myelin.maps$label
myelin.map.PC1$PC1 <- myelin.map.PC1$PC1*-1 #sign of loadings is arbitrary; flip for positive correlation
myelin.map.PC1 %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = PC1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-1.2, 3.5), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/PC1_myelinmaps_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.9, height = 2.7)
myelin.maps <- merge(myelin.maps, myelin.map.PC1, sort = F)

write.csv(myelin.maps, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/myelin_maps.csv", quote = F, row.names = F)

R1-PC1 Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$PC1, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$PC1
## S = 2304070, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6355542

R1-PC1 Correlation Plot

myelin.maps %>%
ggplot(., aes(x = mean.R1, y = PC1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nR1", y="Myelin Map PC1\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/PC1_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.65, height = 1.6)

R-squared

myelin.maps.regression <- lm(scale(myelin.maps$PC1) ~ scale(myelin.maps$mean.R1))
summary(myelin.maps.regression)$r.squared
## [1] 0.4385517

R1 Indexes Cortical Myelin Content Throughout the Cortical Ribbon

Depth-dependent R1 profiles

Developmental study sample

#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
  meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
  meanmap$orig_parcelname <- row.names(meanmap)
  meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
  return(meanmap)
}

regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
  regional_means(depth)
})

regional.myelin.alldepths <-  do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- as.numeric(substr(row.names(regional.myelin.alldepths), 1, 1)) #assign depth variable

regional.myelin.alldepths <-  regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

Across-participant average profiles

area4.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "4_ROI")
area46.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "46_ROI")
areaa24.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "a24_ROI")

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = region, color = region), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.55)

Individual profiles

#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <-  do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- as.numeric(substr(row.names(myelin.glasser.7T.alldepths), 1, 1)) #assign depth variable

#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")

#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")

#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)

#Plot for the oldest 3 participants 
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)
ggplot() +
  geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_individual.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

MICA-PNI high resolution sample

#R1 values in all regions at all depths for all participants
myelin.glasser.micapni.alldepths <- do.call(rbind, myelin.glasser.micapni)
regional.myelin.alldepths.micapni <- myelin.glasser.micapni.alldepths %>% group_by(depth) %>% summarise(across(1:360, \(x) mean(x, na.rm = TRUE))) #average depth R1 per region 
regional.myelin.alldepths.micapni <- pivot_longer(regional.myelin.alldepths.micapni, cols = 2:361, names_to = "label", values_to = "mean.R1")
regional.myelin.alldepths.micapni <- left_join(regional.myelin.alldepths.micapni, glasser.regions, by = "label")

regional.myelin.alldepths.micapni <-  regional.myelin.alldepths.micapni[!(regional.myelin.alldepths.micapni$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

Across-participant average profiles

area4.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "4_ROI")
area46.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "46_ROI")
areaa24.myelin.alldepths <- regional.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(mean.R1 = mean(.$mean.R1)) %>% unnest(cols = c("mean.R1")) %>% mutate(region = "a24_ROI")

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = as.numeric(depth), group = region, color= region), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  scale_x_continuous(breaks = c(0.5, 0.55, 0.6, 0.65)) +
  scale_y_reverse(breaks = c(2, 4, 6, 8, 10, 12,  14), limits = c(14, 2)) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_smooth()`).

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Area_profiles_group_micapni.pdf", device = "pdf", dpi = 500, width = 2.3, height = 1.55)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_smooth()`).

R1 profile coefficient of variation

Developmental study sample

R1.ribbon.cov <- regional.myelin.alldepths %>% group_by(label) %>% do(profile.cov.7T = EnvStats::cv(.$mean.R1)) %>% unnest(cols = c("profile.cov.7T"))

MICA-PNI high resolution sample

R1.ribbon.cov.micapni <- regional.myelin.alldepths.micapni %>% group_by(label) %>% do(profile.cov.micapni = EnvStats::cv(.$mean.R1)) %>% unnest(cols = c("profile.cov.micapni"))

coefficient of variation correlation

R1.ribbon.cov <- left_join(R1.ribbon.cov, R1.ribbon.cov.micapni, by = "label")
write.csv(R1.ribbon.cov, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/R1_cov.csv", quote = F, row.names = F)
cor.test(R1.ribbon.cov$profile.cov.7T, R1.ribbon.cov$profile.cov.micapni, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  R1.ribbon.cov$profile.cov.7T and R1.ribbon.cov$profile.cov.micapni
## S = 2126846, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6635866

R1 profile skewness

Developmental study sample

R1.ribbon.skewness <- regional.myelin.alldepths %>% group_by(label) %>% do(profile.skewness.7T = skewness(.$mean.R1)) %>% unnest(cols = c("profile.skewness.7T"))

MICA-PNI high resolution sample

R1.ribbon.skewness.micapni <- regional.myelin.alldepths.micapni %>% group_by(label) %>% do(profile.skewness.micapni = skewness(.$mean.R1)) %>% unnest(cols = c("profile.skewness.micapni"))

Skewness correlation

R1.ribbon.skewness <- left_join(R1.ribbon.skewness, R1.ribbon.skewness.micapni, by = "label")
write.csv(R1.ribbon.skewness, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/R1_skewness.csv", quote = F, row.names = F)
cor.test(R1.ribbon.skewness$profile.skewness.7T, R1.ribbon.skewness$profile.skewness.micapni, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  R1.ribbon.skewness$profile.skewness.7T and R1.ribbon.skewness$profile.skewness.micapni
## S = 2230618, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6471725

Skewness correlation plot

R1.ribbon.skewness %>%
ggplot(., aes(x = profile.skewness.7T, y = profile.skewness.micapni)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.5, 1.5), breaks = c(-0.5, 0, 0.5, 1, 1.5)) +
  scale_y_continuous(limits = c(-0.5, 1.5), breaks = c(-0.5, 0, 0.5, 1, 1.5)) +
  labs(x="\nSkewness: Developmental", y="Skewness: High-resolution\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure1/Skewness_correlationplot.pdf", device = "pdf", width = 2.68, height = 1.6)

R1-indexed Myelin is Dissociable in Superficial, Middle, and Deep Cortex

Cortical thickness of the frontal lobe

source("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/surface_metrics/surface_measures/extract_surfacestats.R") 
CT.glasser <- extract_surfacestats(myatlas = "glasser", mymeasure = "ThickAvg") #cortical thickness in glasser regions for each scan in the final study sample, generated by fs-tabulate pipeline
CT.glasser.frontal <- CT.glasser[, names(CT.glasser) %in% glasser.frontal$orig_parcelname]
CT.glasser.frontal <- CT.glasser.frontal[, !(names(CT.glasser.frontal) %in% glasser.snr.exclude$orig_parcelname)]
CT.glasser.frontal$subject_id <- CT.glasser$subject_id
CT.glasser.frontal$session_id <- CT.glasser$session_id
CT.glasser.frontal$age <- CT.glasser$age

Across-participant average frontal thickness

average.CT.frontal <- CT.glasser.frontal %>% select(contains("ROI")) %>% colMeans() %>% cbind() %>% as.data.frame() %>% purrr::set_names("CT") %>% mutate(orig_parcelname = row.names(.)) %>% merge(., glasser.regions)
mean(average.CT.frontal$CT)
## [1] 2.613618
sd(average.CT.frontal$CT)
## [1] 0.2025069
ggplot(average.CT.frontal, aes(x = CT)) +
  geom_density(linewidth = 0.5) +
  theme_classic() +
  labs(x="\nCortical Thickness", y="Density\n") +
  theme(
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

Frontal thickness by age bin

region_cols <- setdiff(names(CT.glasser.frontal), c("subject_id", "session_id", "age"))
age.breaks <- seq(10, 32.5, length.out = 6)

average.CT.frontal.agebinned <- CT.glasser.frontal %>%
  mutate(age_bin = cut(age, breaks = age.breaks, include.lowest = TRUE)) %>% 
  group_by(age_bin) %>%
  summarise(
    mean.age = mean(age, na.rm = TRUE),
    n = dplyr::n(),
    across(all_of(region_cols), ~ mean(.x, na.rm = TRUE)),
    .groups = "drop"
)

average.CT.frontal.agebinned <- average.CT.frontal.agebinned %>% pivot_longer(-c(age_bin, mean.age, n), names_to = "region", values_to = "CT")
bin1 <- average.CT.frontal.agebinned %>% filter(age_bin == "[10,14.5]")
bin2 <- average.CT.frontal.agebinned %>% filter(age_bin == "(14.5,19]")
bin3 <- average.CT.frontal.agebinned %>% filter(age_bin == "(19,23.5]")
bin4 <- average.CT.frontal.agebinned %>% filter(age_bin == "(23.5,28]")
bin5 <- average.CT.frontal.agebinned %>% filter(age_bin == "(28,32.5]")

ggplot() +
  geom_density(data = bin5, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[5] ) +
  geom_density(data = bin4, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[4]) +
  geom_density(data = bin3, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[3]) +
  geom_density(data = bin2, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[2]) +
  geom_density(data = bin1, aes(x = CT), linewidth = 0.2, fill = myelin.colorbar[1]) +
  theme_classic() +
  labs(x="\nCortical thickness (mm)", y="Density\n") +
  theme(
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Frontal_CT_agebins.pdf", device = "pdf", width = 2.1, height = 1.8)

Independent voxel compartments in the frontal lobe

ribbon_voxelcounts <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/ribbon_voxelcounts/ribbon_FS_voxelcounts.RDS") #voxel count distributions in the frontal lobe for each longitudinal study session
ribbon_voxelcounts <- left_join(participants, ribbon_voxelcounts, by = c("subject_id", "session_id")) #now for only our final study sample
#Identify the most frequent number of compartments in the cortical ribbon for all study participants
ribbon.compartmentnum.mode <- ribbon_voxelcounts %>% group_by(subject_id, session_id) %>% do(compartment.number = .$voxelcount[which.max(.$count.percents)]) %>% unnest(cols = c("compartment.number"))

table(ribbon.compartmentnum.mode$compartment.number)
## 
##   3 
## 215

Across-participant average voxel compartment number percentages

#Calculate the percent of vertices that show each number of unique voxel compartments
average.counts.frontal <-  ribbon_voxelcounts %>% group_by(voxelcount) %>% do(mean.percent = mean(.$count.percents)) %>% unnest(cols = c("mean.percent"))
kable(average.counts.frontal)
voxelcount mean.percent
1 1.9434674
2 25.3357269
3 45.2212544
4 23.6069641
5 3.7553686
6 0.1370022
7 0.0011926
#98% of vertices have 2 or more unique voxels
average.counts.frontal %>% filter(voxelcount > 1) %>% colSums(.)
##   voxelcount mean.percent 
##     27.00000     98.05751
#73% of vertices have 3 or more unique voxels
average.counts.frontal %>% filter(voxelcount > 2) %>% colSums(.)
##   voxelcount mean.percent 
##     25.00000     72.72178
#Cumulative percentages plot
average.counts.frontal <- average.counts.frontal %>%
  arrange(desc(voxelcount)) %>%
  mutate(cumulative_percentage = cumsum(mean.percent)) %>%
  arrange(voxelcount)

average.counts.frontal %>% filter(voxelcount < 6) %>% #exclude .1% of vertices with 6+ 
ggplot(., aes(x = cumulative_percentage, y = voxelcount, xend = 0)) +
  geom_segment(color = myelin.colorbar[4]) +
  geom_point(size = 1) +
  theme_classic() +
  labs(x="\nCumulative percentage", y="Unique compartments\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Cumulativepercentage_compartments.pdf", device = "pdf", width = 1.7, height = 1.8)
average.counts.frontal %>% filter(voxelcount < 6) %>% #exclude .1% of vertices with 6+ 
ggplot(., aes(x = voxelcount, y = mean.percent)) +
  geom_col(fill = myelin.colorbar[4]) +
  theme_classic() +
  labs(x="\nUnique Compartments", y="Percent of Frontal Vertices\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

Voxel compartment percentages by age bin

ribbon_voxelcounts.agebinned <- ribbon_voxelcounts %>% filter(voxelcount < 6) %>% mutate(age_bin = cut(age, breaks = age.breaks, include.lowest = TRUE)) %>% 
  group_by(age_bin, voxelcount) %>% 
  do(mean.age = mean(.$age), mean.percent = mean(.$count.percents)) %>% unnest(cols = c("mean.age", "mean.percent"))
bin1 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "[10,14.5]")
bin2 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(14.5,19]")
bin3 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(19,23.5]")
bin4 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(23.5,28]")
bin5 <- ribbon_voxelcounts.agebinned %>% filter(age_bin == "(28,32.5]")

ggplot() +
  geom_col(data = bin5, aes(x = voxelcount, y = mean.percent), width = .8, fill = myelin.colorbar[5], se = F ) +
  geom_col(data = bin4, aes(x = voxelcount, y = mean.percent), width = 0.65, fill = alpha(myelin.colorbar[4], 1), se = F) +
  geom_col(data = bin3, aes(x = voxelcount, y = mean.percent), width = 0.5, fill = alpha(myelin.colorbar[3], 1), se = F) +
  geom_col(data = bin2, aes(x = voxelcount, y = mean.percent), width = 0.35, fill = alpha(myelin.colorbar[2], 1), se = F) +
  geom_col(data = bin1, aes(x = voxelcount, y = mean.percent), width = 0.2, fill = alpha(myelin.colorbar[1], 1), se = F) +
  theme_classic() +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5)) +
  scale_y_continuous(breaks = c(0, 10, 20, 30, 40, 50)) +
  labs(x="\nUnique compartments", y="Percent of frontal vertices\n") +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Frontal_compartments_agebins.pdf", device = "pdf", width = 2.1, height = 1.8)

Demarcating superficial, middle, and deep compartments in frontal lobe based on R1 depth derivative

MICA-PNI high resolution sample

#Average depth-dependent R1 profile for the frontal lobe in high resolution (.5 mm) data
frontal.myelin.alldepths.micapni <- regional.myelin.alldepths.micapni[regional.myelin.alldepths.micapni$orig_parcelname %in% glasser.frontal$orig_parcelname,]
depthR1.micapni <- frontal.myelin.alldepths.micapni %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1"))

Frontal lobe R1 depth profile

ggplot(data = depthR1.micapni, aes(x = depth.R1, y = depth)) +
  geom_point() +
  geom_smooth(linewidth = .8, se = F, method = "gam", formula =y ~ s(x, k = 7), color = myelin.colorbar[4]) +
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

Second derivative plot for demarcating compartment boundaries

#Fit a gam to the depth-by-R1 data and calculate the second derivative (rate of change in the rate of change to identify transition zones)
depthR1.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.micapni) 
depthR1.derivatives.micapni <- derivatives(depthR1.model.micapni, order = 2, eps = .1, interval = "simultaneous") #second derivative

ggplot() +
  geom_path(data = depthR1.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.3) +
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.003, 0.006)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/R1_secondderivative_micapni.pdf", device = "pdf", width = 1.6, height = 1.8)
#Repeat the above for our exemplar regions to see if the same transition zones emerge
depthR1.4.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "4_ROI")
depthR1.4.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.4.micapni) 
depthR1.4.derivatives.micapni <- derivatives(depthR1.4.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

depthR1.46.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "46_ROI")
depthR1.46.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.46.micapni) 
depthR1.46.derivatives.micapni <- derivatives(depthR1.46.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

depthR1.a24.micapni <- frontal.myelin.alldepths.micapni %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI") %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1")) %>% mutate(region = "a24_ROI")
depthR1.a24.model.micapni <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.a24.micapni) 
depthR1.a24.derivatives.micapni <- derivatives(depthR1.a24.model.micapni, order = 2, eps = .1, interval = "simultaneous") 

ggplot() +
  geom_path(data = depthR1.a24.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[2]) +
  geom_path(data = depthR1.46.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[3]) +
  geom_path(data = depthR1.4.derivatives.micapni, aes(x = abs(derivative), y = data), linewidth = 0.2, color = myelin.colorbar[5]) +
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.003, 0.006)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2_supplementary/R1_secondderivative_micapni_exemplarregions.pdf", device = "pdf", width = 2.2, height = 2.4)

Change point analysis to validate compartment boundaries

#Identify change points in the relationship between 2 variables. Segmented creates broken-line regression models where relationship between x and y is piece-wise linear and broken up via break points. It identifies where there is a break in the relationship between x and y 
depthR1.lm.micapni = lm(depth.R1 ~ depth, data = depthR1.micapni)
break.points <- segmented::segmented(depthR1.lm.micapni, seg.Z = ~depth, npsi = 2) #identify to depth break points via change point analysis
break.points$indexU
## $depth
##  U1.depth  U2.depth 
##  3.230989 11.400669

Developmental study sample

#Average depth-dependent R1 profile for the frontal lobe in 1 mm data
frontal.myelin.alldepths <- regional.myelin.alldepths[regional.myelin.alldepths$orig_parcelname %in% glasser.frontal$orig_parcelname,]
depthR1.7T <- frontal.myelin.alldepths %>% group_by(depth) %>% do(depth.R1 = mean(.$mean.R1)) %>% unnest(cols = c("depth.R1"))
ggplot(data = depthR1.7T, aes(x = depth.R1, y = depth)) +
  geom_point() +
  geom_smooth(linewidth = .8, method = "gam", formula = y ~ s(x, k = 6), color = myelin.colorbar[4]) + 
  theme_classic() +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

#Fit a gam to the depth-by-R1 data and calculate the second derivative
depthR1.model <- gam(depth.R1 ~ s(depth, k = 7), method = "REML", data = depthR1.7T)
depthR1.derivatives <- derivatives(depthR1.model, order = 2, eps = .1, interval = "simultaneous") #second derivative

ggplot(depthR1.derivatives, aes(x = abs(derivative), y = data)) +
  geom_path(linewidth = 0.3) + 
  theme_classic() +
  scale_x_continuous(breaks = c(0, 0.005, 0.01)) +
  scale_y_reverse(breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  xlab("\nSecond derivative") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2_supplementary/R1_secondderivative_7T.pdf", device = "pdf", width = 2.2, height = 2.4)

Regional R1 Distributions in Superficial, Middle, Deep

myelin.smd.7T <- do.call(rbind, myelin.compartments.7T) #R1 in superficial, middle, deep for all participants and regions

myelin.smd.7T <- myelin.smd.7T[ , 
  names(myelin.smd.7T) %in% glasser.frontal$orig_parcelname] #just frontal
myelin.smd.7T <- myelin.smd.7T[ , 
  !(names(myelin.smd.7T) %in% glasser.snr.exclude$orig_parcelname)] #good SNR

myelin.smd.7T$depth <- str_remove(row.names(myelin.smd.7T), "\\..*")

cols_to_pivot <- names(myelin.smd.7T)[grepl("ROI", names(myelin.smd.7T))]
myelin.smd.7T <- myelin.smd.7T %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1") #long format
myelin.smd.7T <- left_join(myelin.smd.7T, glasser.regions, by = "orig_parcelname")

myelin.smd.7T <- myelin.smd.7T %>% group_by(depth, label) %>% do(mean.depth.R1 = mean(.$R1)) %>% unnest(cols = c("mean.depth.R1")) #R1 in superficial, middle, deep in frontal regions averaged across participants
myelin.smd.7T %>% 
  ggplot(aes(x = mean.depth.R1, y = depth)) +
  geom_violin(aes(fill = depth), color = "white", bw = .006) +
  stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
  theme_classic() +
  scale_fill_manual(values = c("#004A38FF", "#809dc4", "#A29DC4")) +
  scale_x_continuous(limits = c(0.497, 0.65), expand = c(0,0)) +
  scale_y_discrete(expand = c(0, 0)) +
  xlab("\nR1") +
  ylab("Cortical depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/Depthplot_regionalR1_frontal.pdf", device = "pdf", dpi = 500, width = 2.2, height = 1.7)
for(this.depth in c(unique(myelin.smd.7T$depth))){
  
  depth.data <- myelin.smd.7T %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  depth.plot <- ggplot() + 
  geom_brain(data = depth.data, atlas = glasser,  mapping = aes(fill = mean.depth.R1), colour=I("gray10"), size=I(.05)) + theme_classic() + theme(legend.position = "none") + 
   paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.48, 0.6), oob = squish, na.value = "white") + theme(legend.position = "none") +
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray10", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

  print(depth.plot)
  ggsave(filename = sprintf("/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure2/%s_corticalmap_frontal.pdf", this.depth), device = "pdf", dpi = 500, width = 4.1, height = 1.5)
}